#include "Riostream.h"
#include <iostream>
#include <iomanip>
#include <string>

#include "TFile.h"
#include "TTree.h"
#include "TChain.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TCanvas.h"
#include "TPad.h"
#include "TStyle.h"
#include "TROOT.h"
#include "TMath.h"
#include "TRandom.h"
#include "TFormula.h"
#include "TSystem.h"
#include "TEnv.h"

#include "CommandLineInterface.hh"
#include "Barrel.hh"
#include "Annular.hh"
#include "Germanium.hh"
#include "SRIMloader.hh"

using namespace TMath;
using namespace std;

ClassImp(Barrel);
ClassImp(Annular);
ClassImp(Germanium);

SRIM srim; 

int main(int argc, char* argv[]){
  vector<char*> InputFiles;
  char* OutputFile = NULL;
  char* MBFile = NULL;
  bool verbose = false;
  int all = 1;
  int tiltcorr = 0;
  
  Float_t z,yprime,offset,r,alpha,theta,phi, thetadet;
  Float_t Rafter, Efor, Ecor, dEcor, Eex, Ein;
  Float_t gamma, beta, EGamDS;
  Float_t m[20], MBTheta[8][3][6];
  Int_t par;

  // Masses
  m[1] = 1.00782503207*931.494e3;
  m[2] = 2.01410177785*931.494e3;
  m[3] = 3.01604927767*931.494e3;
  m[4] = 4.00260325415*931.494e3;
  m[10] = 10.013533818*931.494e3;
  m[11] = 11.021657749*931.494e3;
  m[12] = 12.026920737*931.494e3;


  // Beam energy
  Ein = 2.65e3*m[11]/931.434e3;


  CommandLineInterface* interface = new CommandLineInterface();

  interface->Add("-i", "inputfiles", &InputFiles);
  interface->Add("-o", "outputfile", &OutputFile);
  interface->Add("-M", "miniball angle", &MBFile);
  interface->Add("-a", "all histos? 0 only Erear vs Pos for each strip 2 no gammas and no cd", &all);
  interface->Add("-t", "1 tiltcorr histos, 0 no tiltcorr", &tiltcorr);
  interface->Add("-v", "verbose", &verbose);
  
  interface->CheckFlags(argc, argv);

  if(InputFiles.size() == 0 || OutputFile == NULL){
    cerr<<"You have to provide at least one input file and the output file!"<<endl;
    exit(1);
  }

  cout<<"input file(s):"<<endl;
  for(unsigned int i=0; i<InputFiles.size(); i++){
    cout<<InputFiles[i]<<endl;
  }
      
  cout<<"output file: "<<OutputFile<< endl;

  // Generating SRIM graphs for energyloss calculations
  srim.Creategraphs();
//  printf("R = %f, E = %f \n",srim.RMylar[4]->Eval(3.0),srim.EMylar[4]->Eval(11.0));
  Rafter = srim.RT[11]->Eval(Ein);
  Ein = srim.ET[11]->Eval(Rafter - 5);
  gamma = (Ein + m[11])/m[11];
  beta = sqrt(1 - 1/(gamma*gamma));

	TEnv *Angels = new TEnv(MBFile);
	for (int i=0; i<8; i++){
		for (int j=0; j<3; j++){
			for (int k=0; k<6; k++){
				MBTheta[i][j][k] = Angels->GetValue(Form("Theta.%d.%d.%d",i,j,k),0.0);
				printf("theta = %f \n",MBTheta[i][j][k]);
			}
		}
	}

  TChain* tr;
  tr = new TChain("caltr");
  for(unsigned int i=0; i<InputFiles.size(); i++){
    tr->Add(InputFiles[i]);
  }

  if(tr == NULL){
    cout << "could not find tree caltr in file " << endl;
    for(unsigned int i=0; i<InputFiles.size(); i++){
      cout<<InputFiles[i]<<endl;
    }
    return 3;
  }

  vector<Barrel> *FBarrel = new vector<Barrel>;
  vector<Barrel> *BBarrel = new vector<Barrel>;
  vector<Annular> *FCD = new vector<Annular>;
  vector<Germanium> *Miniball = new vector<Germanium>;
  long long EbisTime;
  long long T1Time;
  long long SuperCycleTime;

  tr->SetBranchAddress("ForwardBarrel",&FBarrel);
  tr->SetBranchAddress("BackwardBarrel",&BBarrel);
  tr->SetBranchAddress("ForwardCD",&FCD);
  tr->SetBranchAddress("Miniball",&Miniball);
  tr->SetBranchAddress("EbisTime",&EbisTime);
  tr->SetBranchAddress("T1Time",&T1Time);
  tr->SetBranchAddress("SuperCycleTime",&SuperCycleTime);

  TFile* outfile = new TFile(OutputFile,"recreate");

  if(outfile->IsZombie()){
    return 4;
  }
  
	// Generating cut functions
	TFormula* back1 = new TFormula("back1","3200 - 0.5*x");
	TFormula* back2 = new TFormula("back2","4000 - 0.35*x");
	TFormula* for1 = new TFormula("for1","10e6/(x + 3.5e3)");
	TFormula* for2 = new TFormula("for2","21e6/(x + 4.5e3)");
	TFormula* for3 = new TFormula("for3","33e6/(x + 5.e3)");
	TFormula* for4 = new TFormula("for4","52e6/(x + 7e3)");
	TFormula* CD1 = new TFormula("CD1","27.8 - .05*x");
	TFormula* CD2 = new TFormula("CD2","30. - .05*x");

	
	// Prepare for vectorial histograms
	// Barrel
	// Position
	TH1F* ForwardBarrel_Pos[4][16];
	TH2F* ForwardBarrel_PosvsStrip[4];
	TH1F* BackwardBarrel_Pos[4][16];
	TH2F* BackwardBarrel_PosvsStrip[4];
	// dE
	TH1F* ForwardBarrel_dE[4][16];
	TH1F* ForwardBarrel_dEtot[4];
	TH2F* ForwardBarrel_dEvsStrip[4];
	TH2F* ForwardBarrel_dEvsDetangle[4];
	TH1F* ForwardBarrel_dEangle[4][20];
	TH1F* ForwardBarrel_dEforPos[4][16][16];
	TH2F* ForwardBarrel_dEvsPos[4][16];
	TH1F* BackwardBarrel_dE[4][16];
	TH1F* BackwardBarrel_dEtot[4];
	TH2F* BackwardBarrel_dEvsStrip[4];
	TH1F* BackwardBarrel_dEforPos[4][16];
	TH2F* BackwardBarrel_dEvsPos[4];

	// Efor
	TH1F* ForwardBarrel_Efor[4];
	TH2F* ForwardBarrel_EforvsDetangle[4];
	TH2F* ForwardBarrel_EforvsPos[4][16];

	// E backdetector
	TH2F *ForwardBarrel_dEvsEtot[4];
	TH2F *ForwardBarrel_dEvsE[4][16];
	TH2F* ForwardBarrel_dEvsEangle[4][20];
	TH2F *BackwardBarrel_dEvsEtot[4];
	TH2F *BackwardBarrel_dEvsE[4][16];
	TH2F* BackwardBarrel_dEvsEangle[4][20];
	TH2F* ForwardBarrel_dEvsEcor[4];
	TH2F* ForwardBarrel_dEvsEcord[4];

	//cd
	// dE
	TH1F* CD_dEtot[4];
	TH2F* CD_dEvsStrip[4];
	TH2F* CD_dEvsRing[4];
	TH2F* CD_StripvsRing[4];
	TH2F* CD_EStripvsERing[4];
	TH1F* CD_dEforStrip[4][16];
	TH1F* CD_dEforRing[4][16];
	
	// E backdetector
	TH2F *CD_dEvsEtot[4];
	TH2F *CD_dEvsEStrip[4][16];
	TH2F *CD_dEvsERing[4][16];
	
	// Angle
	TH2F *CD_dEvsTheta[4];
	TH1F *Be11_Theta[4];

	// Gammas
	// Energy
	TH1F* Gamma_ECrystals[8][3];
	TH1F* Gamma_ESegments[8][3][6];
	TH1F* Gamma_Protons_ESegments[8][3][6];
	TH1F* Gamma_Deuterons_ESegments[8][3][6];
	

	// Particles
	TH2F* Protons_EtotvsLabangle[4];
	TH2F* Deuterons_EtotvsLabangle[4];
	TH2F* Tritons_EtotvsLabangle[4];
	TH1F* Protons_Eex[4];
	TH1F* Deuterons_Eex[4];
	TH1F* Tritons_Eex[4];


	// No punchthrough
	TH2F* ForwardBarrel_dEvsLabangle[4];
	TH2F* BackwardBarrel_dEvsLabangle[4];

	// Initiate histograms
	
	// Single histograms
	// Angle histograms
	TH1F* ForwardBarrel_Theta = new TH1F("ForwardBarrel_Theta","Theta in forward barrel",100,0,TMath::Pi()/2.);
	TH1F* ForwardBarrel_Phi = new TH1F("ForwardBarrel_Phi","Phi in forward barrel",100,0.,2*TMath::Pi());
	TH1F* ForwardBarrel_Theta_det = new TH1F("ForwardBarrel_Theta_det","detectorangles for forward barrel",100,0,TMath::Pi()/2.);
	TH1F* BackwardBarrel_Theta = new TH1F("BackwardBarrel_Theta","Theta in backward barrel",100,0,TMath::Pi()/2.);
	TH1F* BackwardBarrel_Phi = new TH1F("BackwardBarrel_Phi","Phi in backward barrel",100,0.,2*TMath::Pi());
	TH1F* CD_Theta = new TH1F("CD_Theta","Theta in CD",100,0,TMath::Pi()/2.);
	TH1F* CD_Phi = new TH1F("CD_Phi","Phi in CD",100,0.,2*TMath::Pi());

	// E backdetector
	TH2F *ForwardBarrel_dEvsEtotal = new TH2F("ForwardBarrel_dEvsEtotal","dE vs Eback for fordet",500,-1e3,1.5e4,500,0,1e4);
	TH2F *BackwardBarrel_dEvsEtotal = new TH2F("BackwardBarrel_dEvsEtotal","dE vs Eback for backdet",500,-1e3,1e4,500,0,1e4);
	TH2F *ForwardBarrel_dEvsEcortot = new TH2F("ForwardBarrel_dEvsEcortot","dE vs Eback corrected for fordet",500,0,15.,500,0,10.);
	TH2F *BackwardBarrel_dEvsEcortot = new TH2F("BackwardBarrel_dEvsEcortot","dE vs Eback corrected for backdet",500,0,1e4,500,0,1e4);
	TH2F *ForwardBarrel_dEvsEcorp = new TH2F("ForwardBarrel_dEvsEcorp","dE vs Eback corrected protons for fordet",500,0,1.5e4,500,0,1e4);
	for (Int_t i=0;i<4;i++){ForwardBarrel_dEvsEcord[i] = new TH2F(Form("ForwardBarrel_dEvsEcord_%d",i),Form("dE vs Eback corrected deuterons for fordet %f",i),500,0,1.5e4,500,0,1e4);}
	TH2F *ForwardBarrel_dEvsEcort = new TH2F("ForwardBarrel_dEvsEcort","dE vs Eback corrected tritons for fordet",500,0,1.5e4,500,0,1e4);
	TH2F *BackwardBarrel_dEvsEcorp = new TH2F("BackwardBarrel_dEvsEcorp","dE vs Eback corrected protons for backdet",500,0,1e4,500,0,1e4);

	// Particle histograms
	TH2F* Protons_EtotvsLabangletot = new TH2F("Protons_EtotvsLabangletot","E vs labangle for protons",180,0,180,500,0,1.6e4);
	TH2F* Deuterons_EtotvsLabangletot = new TH2F("Deuterons_EtotvsLabangletot","E vs labangle for deuterons",180,0,90,500,0,1.6e4);
	TH2F* Tritons_EtotvsLabangletot = new TH2F("Tritons_EtotvsLabangletot","E vs labangle for tritons",180,0,90,500,0,2e4);
	TH1F* Protons_Eextot = new TH1F("Protons_Eextot","Exitation energy for protons in fordet",200,-2.,4.);
	TH1F* Deuterons_Eextot = new TH1F("Deuterons_Eextot","Exitation energy for deuterons in fordet",200,-2,3.5);
	TH1F* Tritons_Eextot = new TH1F("Tritons_Eextot","Exitation energy for tritons in fordet",200,-2,10.);
	TH2F* CD_EtotvsLabangle = new TH2F("CD_EtotvsLabangle","E vs labangle in CD",360,0,30,2e3,0,3e4);

	// No punchthrough
	TH2F* ForwardBarrel_dEvsLabangletot = new TH2F("ForwardBarrel_dEvsLabangletot","E vs labangle in for de detector",180,0,90,500,0,15.);
	TH2F* ForwardBarrel_dEvsLabangle_gammacut = new TH2F("ForwardBarrel_dEvsLabangle_gammacut","E vs labangle in for de detector gated on neutrons",180,0,90,500,0,15.);

	// CD
	// Angle
	TH2F* CD_dEvsThetatot = new TH2F("CD_dEvsThetatot","dE vs Theta for CD",90,0,35,1e4,0,35.);

	// Gamma
	TH1F* Gamma_Etot = new TH1F("Gamma_Etot","Gamma energies",1000,0,10e3);
	TH1F* Gamma_EProton = new TH1F("Gamma_EProton","Gamma energies",400,0,3e3);
	TH1F* Gamma_EDeuteron = new TH1F("Gamma_EDeuteron","Gamma energies",400,0,2e3);
	TH1F* Gamma_ETriton = new TH1F("Gamma_ETriton","Gamma energies",400,0,3.5e3);
	TH1F* Gamma_ENPT = new TH1F("Gamma_ENPT","Gamma energies",400,0,4e3);
	TH1F* Gamma_EProtonDS = new TH1F("Gamma_EProtonDS","Gamma energies",400,0,3e3);
	TH1F* Gamma_EDeuteronDS = new TH1F("Gamma_EDeuteronDS","Gamma energies",400,0,2e3);
	TH1F* Gamma_ETritonDS = new TH1F("Gamma_ETritonDS","Gamma energies",400,0,3.5e3);
	TH1F* Gamma_ENPTDS = new TH1F("Gamma_ENPTDS","Gamma energies",400,0,4e3);
	TH1F* Gamma_EDoppler = new TH1F("Gamma_EDoppler","Gamma energies",300,0,4e3);
	TH2F* Gamma_EexvsEGammaProton = new TH2F("Gamma_EexvsEGammaProton","Eex vs Egamma for 12Be",100,0,3.5,100,0,3.e3);
	TH2F* Gamma_EexvsEGammaDeuteron = new TH2F("Gamma_EexvsEGammaDeuteron","Eex vs Egamma for 11Be",100,0,1.5,100,0,1.5e3);
	TH2F* Gamma_EexvsEGammaTriton = new TH2F("Gamma_EexvsEGammaTriton","Eex vs Egamma for 10Be",100,0,7.5,100,0,3.e3);

	// Time
	TH2F* Gamma_EvsdtProton = new TH2F("Gamma_EvsdtProton","Gamma energy vs timediff between gamma and proton",60,-400,1200,1e3,0,4.);
	TH2F* Gamma_EvsdtDeuteron = new TH2F("Gamma_EvsdtDeuteron","Gamma energy vs timediff between gamma and deuteron",60,-400,1200,1e3,0,4.);
	TH2F* Gamma_EvsdtTriton = new TH2F("Gamma_EvsdtTriton","Gamma energy vs timediff between gamma and triton",60,-400,1200,1e3,0,4.);
	TH2F* Gamma_EvsdtNPT = new TH2F("Gamma_EvsdtNPT","Gamma energy vs timediff between gamma and NPT particles",60,-400,1200,1e3,0,4.);
	TH2F* Gamma_EvsdtBack = new TH2F("Gamma_EvsdtBack","Gamma energy vs timediff between gamma and particles in back",60,-400,1200,1e3,0,4.);
	TH2F* Gamma_EvsdtCD = new TH2F("Gamma_EvsdtCD","Gamma energy vs timediff between gamma and particles in CD",60,-400,1200,1e3,0,4.);
	TH2F* Gamma_Evsdttot = new TH2F("Gamma_Evsdttot","Gamma energy vs timediff between gamma and particle",240,-3200,3200,1e3,0,4.);
	TH1F* Gamma_dtGe = new TH1F("Gamma_dtGe","Decay time of excited Ge state",60,-400,1200);
	TH1F* Gamma_Ecut = new TH1F("Gamma_Ecut","E gamma cut for delayed decay",1e3,2.e3,2.5e3);

	// Pad dependent histograms
	for (Int_t i=0;i<4;i++){
		// Position
		ForwardBarrel_PosvsStrip[i] = new TH2F(Form("ForwardBarrel_PosvsStrip_%d",i),Form("Position vs stripnumber for fordet %d",i),16,0,16,200,0,1.1);
		BackwardBarrel_PosvsStrip[i] = new TH2F(Form("BackwardBarrel_PosvsStrip_%d",i),Form("Position vs stripnumber for backdet %d",i),16,0,16,200,0,1.1);

		// dE
		ForwardBarrel_dEtot[i] = new TH1F(Form("ForwardBarrel_dEtot_%d",i),Form("Total dE for fordet %d",i),500,0,5.e3);
		ForwardBarrel_dEvsStrip[i] = new TH2F(Form("ForwardBarrel_dEvsStrip_%d",i),Form("dE vs stripnumber for fordet %d",i),16,0,16,300,0,5.e3);
		ForwardBarrel_dEvsDetangle[i] = new TH2F(Form("ForwardBarrel_dEvsDetangle_%d",i),Form("dE vs detector angle for fordet %d",i),100,0,TMath::Pi()/2.,500,0,5.);
		BackwardBarrel_dEtot[i] = new TH1F(Form("BackwardBarrel_dEtot_%d",i),Form("Total dE for backdet %d",i),500,0,1e4);
		BackwardBarrel_dEvsStrip[i] = new TH2F(Form("BackwardBarrel_dEvsStrip_%d",i),Form("dE vs stripnumber for backdet %d",i),16,0,16,300,0,1e4);
		BackwardBarrel_dEvsPos[i] = new TH2F(Form("BackwardBarrel_dEvsPos_%d",i),Form("dE vs position for backdet %d",i),300,0,1.1,300,0,1e4);

		// Efor
		ForwardBarrel_Efor[i] = new TH1F(Form("ForwardBarrel_Efor_%d",i),Form("Total Efor for fordet %d",i),500,4e3,7e3);
		ForwardBarrel_EforvsDetangle[i] = new TH2F(Form("ForwardBarrel_EforvsDetangle_%d",i),Form("Efor vs detector angle for fordet %d",i),500,0,TMath::Pi()/2.,500,0.e3,1.e4);

		// E backdetector
		ForwardBarrel_dEvsEtot[i] = new TH2F(Form("ForwardBarrel_dEvsEtot_%d",i),Form("dE vs Eback for fordet %d",i),500,0.e3,14.e3,300,0,1e4);
		BackwardBarrel_dEvsEtot[i] = new TH2F(Form("BackwardBarrel_dEvsEtot_%d",i),Form("dE vs Eback for backdet %d",i),300,0,14.e3,300,0,7e3);
		ForwardBarrel_dEvsEcor[i] = new TH2F(Form("ForwardBarrel_dEvsEcor_%d",i),Form("dE vs Eback corrected for fordet %d",i),500,0,1.5e4,500,0,1e4);

		// No punchthrough
		ForwardBarrel_dEvsLabangle[i] = new TH2F(Form("ForwardBarrel_dEvsLabangle_%d",i),Form("dE vs labangle for fordet %d",i),180,0,90,500,0,10.e3);
		BackwardBarrel_dEvsLabangle[i] = new TH2F(Form("BackwardBarrel_dEvsLabangle_%d",i),Form("dE vs labangle for backdet %d",i),180,90,180,500,0,10.e3);

		// particles
		Protons_EtotvsLabangle[i] = new TH2F(Form("Protons_EtotvsLabangle_%d",i),Form("E vs labangle for protons in det %d",i),180,0,90,500,0,1.6e4);
		Deuterons_EtotvsLabangle[i] = new TH2F(Form("Deuterons_EtotvsLabangle_%d",i),Form("E vs labangle for deuterons in det %d",i),180,0,90,500,0,1.6e4);
		Tritons_EtotvsLabangle[i] = new TH2F(Form("Tritons_EtotvsLabangle_%d",i),Form("E vs labangle for tritons in det %d",i),180,0,90,500,0,2e4);
		Protons_Eex[i] = new TH1F(Form("Protons_Eex_%d",i),Form("Exitation energy for protons in fordet %d",i),100,-2000,4.e3);
		Deuterons_Eex[i] = new TH1F(Form("Deuterons_Eex_%d",i),Form("Exitation energy for deuterons in fordet %d",i),100,-2000,4.e3);
		Tritons_Eex[i] = new TH1F(Form("Tritons_Eex_%d",i),Form("Exitation energy for tritons in fordet %d",i),100,-1500,1e4);

		// Strip dependent histograms
		for (Int_t j=0;j<16;j++){
			// Position
			ForwardBarrel_Pos[i][j] = new TH1F(Form("ForwardBarrel_Pos_%d_%d",i,j),Form("Position for strip %d in fordet %d",j,i),500,0.,1.1);
			BackwardBarrel_Pos[i][j] = new TH1F(Form("BackwardBarrel_Pos_%d_%d",i,j),Form("Position for strip %d in backdet %d",j,i),500,0.,1.1);
		
			// dE
			ForwardBarrel_dE[i][j] = new TH1F(Form("ForwardBarrel_dE_%d_%d",i,j),Form("dE for strip %d in fordet %d",j,i),500,0,5.e3);
			ForwardBarrel_dEvsPos[i][j] = new TH2F(Form("ForwardBarrel_dEvsPos_%d_%d",i,j),Form("dE vs pos in strip %d for fordet %d",j,i),300,0,1.1,500,0.,5.e3);
			// And position dependent
			for (Int_t k=0;k<16;k++) ForwardBarrel_dEforPos[i][j][k] = new TH1F(Form("ForwardBarrel_dEforPos_%d_%d_%d",i,j,k),Form("dE for pos %d in strip %d in fordet %d",k,j,i),500,0,5.e3);
			
			BackwardBarrel_dE[i][j] = new TH1F(Form("BackwardBarrel_dE_%d_%d",i,j),Form("dE for strip %d in backdet %d",j,i),500,0,1e4);
			BackwardBarrel_dEforPos[i][j] = new TH1F(Form("BackwardBarrel_dEforPos_%d_%d",i,j),Form("dE for pos %d in backdet %d",j,i),500,0,1e4);
		
			// E backdetector
			ForwardBarrel_dEvsE[i][j] = new TH2F(Form("ForwardBarrel_dEvsE_%d_%d",i,j),Form("dE vs Eback for strip %d in fordet %d",j,i),250,0,3.5e3,300,0,1e4);
			BackwardBarrel_dEvsE[i][j] = new TH2F(Form("BackwardBarrel_dEvsE_%d_%d",i,j),Form("dE vs Eback for strip %d in backdet %d",j,i),300,0,1e4,300,0,1e4);
	
			// Efor
			ForwardBarrel_EforvsPos[i][j] = new TH2F(Form("ForwardBarrel_EforvsPos_%d_%d",i,j),Form("Efor vs pos in strip %d for fordet %d",j,i),300,0,1.1,500,0,9.e3);
		}
	
		// Angle depent histograms
		for (Int_t j =0;j<20;j++){
			// dEfile
			ForwardBarrel_dEangle[i][j] = new TH1F(Form("ForwardBarrel_dEangle_%d_%d",i,j),Form("dE around angle %d in fordet %d",j,i),500,0,1.5e3);
			ForwardBarrel_dEvsEangle[i][j] = new TH2F(Form("ForwardBarrel_dEvsEangle_%d_%d",i,j),Form("dE vs E around angle %d in fordet %d",j,i),250,0.,10.,250,0,10.);
			BackwardBarrel_dEvsEangle[i][j] = new TH2F(Form("BackwardBarrel_dEvsEangle_%d_%d",i,j),Form("dE vs E around angle %d in backdet %d",j,i),500,0,1e4,500,0,1e4);
		}
	}
	
	// CD dependent histograms
	for(int i=0;i<4;i++){
		CD_dEtot[i] = new TH1F(Form("CD_dEtot_%d",i), Form("dE for CD %d",i),600,0,2e4);
		CD_dEvsStrip[i] = new TH2F(Form("CD_dEvsStrip_%d",i), Form("deltaE vs Strip Nr of CD %d",i),16,0,16,500,0,2e4);
		CD_dEvsRing[i] = new TH2F(Form("CD_dEvsRing_%d",i), Form("deltaE vs Ring Nr ofCD %d",i),16,0,16,500,0,2e4);
		CD_StripvsRing[i] = new TH2F(Form("CD_StripvsRing_%d",i), Form("Strip vs Ring of CD %d",i),16,0,16,16,0,16);
		CD_EStripvsERing[i] = new TH2F(Form("CD_EStripvsERing_%d",i), Form("Strip vs Ring Energy of CD %d",i),500,0,2e4,500,0,2e4);
		CD_dEvsEtot[i] = new TH2F(Form("CD_dEvsEtot_%d",i), Form("dE vs E for CD %d",i),200,0,100,300,0,2.5e4);

		// Angle
		CD_dEvsTheta[i] = new TH2F(Form("CD_dEvsTheta_%d",i),Form("dE vs Theta for CD %d",i),90,0,35,1e4,0,35.);

		// 11Be
		Be11_Theta[i] = new TH1F(Form("Be11_Theta_%d",i),Form("Theta for 11Be in CD %d",i),90,0,35);

		// Strip/ring dependent histograms
		for(int j=0;j<16;j++){
			CD_dEforStrip[i][j] = new TH1F(Form("CD_dEforStrip_%d_%d",i,j), Form("dE for strip %d in CD %d",j,i), 1e3,0,2e4);
			CD_dEforRing[i][j] = new TH1F(Form("CD_dEforRing_%d_%d",i,j), Form("dE for ring %d in CD %d",j,i), 1e3,0,2e4);
			CD_dEvsEStrip[i][j] = new TH2F(Form("CD_dEvsEStrip_%d_%d",i,j), Form("dE vs E for strip %d in CD %d",j,i),300,0,100,300,0,2e4);
			CD_dEvsERing[i][j] = new TH2F(Form("CD_dEvsERing_%d_%d",i,j), Form("dE vs E for ring %d in CD %d",j,i),300,0,100,300,0,2e4);
		}
	}
	
	// Gammas
	for (Int_t i=0;i<8;i++){
		
		//Crystals
		for (Int_t j=0;j<3;j++){
			Gamma_ECrystals[i][j] = new TH1F(Form("Gamma_ECrystals_%d_%d",i,j),Form("Gamma energies in cluster %d and crystal %d",i,j),1000,0,4e3);

			// Segments
			for (Int_t k=0;k<6;k++){
				Gamma_ESegments[i][j][k] = new TH1F(Form("Gamma_ESegments_%d_%d_%d",i,j,k),Form("Gamma energies in cluster %d,crystal %d and segment %d",i,j,k),1000,800,2.5e3);
				Gamma_Protons_ESegments[i][j][k] = new TH1F(Form("Gamma_Protons_ESegments_%d_%d_%d",i,j,k),Form("Gamma energies in cluster %d,crystal %d and segment %d",i,j,k),250,0,2.5e3);
				Gamma_Deuterons_ESegments[i][j][k] = new TH1F(Form("Gamma_Deuterons_ESegments_%d_%d_%d",i,j,k),Form("Gamma energies in cluster %d,crystal %d and segment %d",i,j,k),250,0,3.5e3);

			}
		}
	}


// END HISTOGRAMS

bool partfound;
  double time[12];
  Double_t nentries = tr->GetEntries();
  Int_t nbytes = 0;
  Int_t status;
  for(int i=0; i<nentries;i++){
    par = 0;
	 status = tr->GetEvent(i);
    if(status == -1){
      cerr<<"Error occured, couldn't read entry "<<i<<" from tree "<<tr->GetName()<<" in file "<<tr->GetFile()->GetName()<<endl;
      return 5;
    }
    else if(status == 0){
      cerr<<"Error occured, entry "<<i<<" in tree "<<tr->GetName()<<" in file "<<tr->GetFile()->GetName()<<" doesn't exist"<<endl;
      return 6;
    }
    nbytes += status;

    partfound =false;
    for(int k=0;k<12;k++) time[k]=0.;

    //forward barrel
    for(unsigned int j=0; j<FBarrel->size(); j++){
		if((*FBarrel)[j].GetMult() ==0) continue;

		if( (*FBarrel)[j].GetID() < 0 || (*FBarrel)[j].GetID() > 3 ){
			cerr<<"Error in entry "<<i<<": "<<j<<". forward barrel detector id is wrong: "<< (*FBarrel)[j].GetID()<<endl;	
			continue;
		}

		if( (*FBarrel)[j].GetStripNr()[0] < 0 || (*FBarrel)[j].GetStripNr()[0] > 15 ){
			cerr<<"Error in entry "<<i<<": "<<(*FBarrel)[j].GetID()<<". forward barrel detector strip nr is wrong: "<<(*FBarrel)[j].GetStripNr()[0]<<endl;	
			continue;
		}
		partfound=true;
		time[(*FBarrel)[j].GetID()] = (*FBarrel)[j].GetTime();
	
		//filling histos
		if(all!=0){
			// Position histograms
			ForwardBarrel_Pos[(*FBarrel)[j].GetID()][(*FBarrel)[j].GetStripNr()[0]]->Fill((*FBarrel)[j].GetStripPos()[0]);
			ForwardBarrel_PosvsStrip[(*FBarrel)[j].GetID()]->Fill((*FBarrel)[j].GetStripNr()[0],(*FBarrel)[j].GetStripPos()[0]);

			// Cylindrical coordinates
			z = 8. + 3.125*((*FBarrel)[j].GetStripNr()[0] + gRandom->Rndm());
			yprime = 29. + 0.*((*FBarrel)[j].GetID()==0) + 2.*((*FBarrel)[j].GetID()==3) - 2.*((*FBarrel)[j].GetID()==1) - 0.*((*FBarrel)[j].GetID()==2);
//			yprime = 29.;
			offset = 0 + .04*((*FBarrel)[j].GetID()==0) + 0.*((*FBarrel)[j].GetID()==1) - .04*((*FBarrel)[j].GetID()==2) - 0.*((*FBarrel)[j].GetID()==3);
//			offset = 0.;
			r = sqrt(yprime*yprime + 50*50*(0.5 + offset - (*FBarrel)[j].GetStripPos()[0])*(0.5 + offset - (*FBarrel)[j].GetStripPos()[0]));
			alpha = atan(((*FBarrel)[j].GetStripPos()[0] - 0.5 - offset)*50/yprime);
			// Spherical coordinates
			theta = atan(r/z) +.005*((*FBarrel)[j].GetID()==0) + 0.*((*FBarrel)[j].GetID()==3) - 0.*((*FBarrel)[j].GetID()==1) - .005*((*FBarrel)[j].GetID()==2);
//			theta = atan(r/z);
			phi = alpha + TMath::Pi()/2.*(*FBarrel)[j].GetID();
			if (phi< 0) phi = phi+2*TMath::Pi();
			// Detector angle
			thetadet = acos(sin(theta)*cos(alpha));
			// angle histograms
			ForwardBarrel_Theta->Fill(theta);
			ForwardBarrel_Phi->Fill(phi);
			ForwardBarrel_Theta_det->Fill(thetadet);

			// dE histograms
			ForwardBarrel_dE[(*FBarrel)[j].GetID()][(*FBarrel)[j].GetStripNr()[0]]->Fill((*FBarrel)[j].GetRear());
			ForwardBarrel_dEvsPos[(*FBarrel)[j].GetID()][(*FBarrel)[j].GetStripNr()[0]]->Fill((*FBarrel)[j].GetStripPos()[0],(*FBarrel)[j].GetRear());
//			ForwardBarrel_dEvsPos[(*FBarrel)[j].GetID()][0]->Fill((*FBarrel)[j].GetStripPos()[0],(*FBarrel)[j].GetRear());
			ForwardBarrel_dEtot[(*FBarrel)[j].GetID()]->Fill((*FBarrel)[j].GetRear());
			ForwardBarrel_dEvsStrip[(*FBarrel)[j].GetID()]->Fill((*FBarrel)[j].GetStripNr()[0],(*FBarrel)[j].GetRear());
			ForwardBarrel_dEvsDetangle[(*FBarrel)[j].GetID()]->Fill(thetadet,(*FBarrel)[j].GetRear()/1.e3);
			for (Int_t M=0;M<16;M++){
				if (((*FBarrel)[j].GetStripPos()[0] > M*.06667)&&((*FBarrel)[j].GetStripPos()[0] < (M+1)*.06667)){
					ForwardBarrel_dEforPos[(*FBarrel)[j].GetID()][(*FBarrel)[j].GetStripNr()[0]][M]->Fill((*FBarrel)[j].GetRear());
				}
			}

			// E back histograms
			ForwardBarrel_dEvsEtotal->Fill((*FBarrel)[j].GetEdet(),(*FBarrel)[j].GetRear());
			ForwardBarrel_dEvsEtot[(*FBarrel)[j].GetID()]->Fill((*FBarrel)[j].GetEdet(),(*FBarrel)[j].GetRear());
			ForwardBarrel_dEvsE[(*FBarrel)[j].GetID()][(*FBarrel)[j].GetStripNr()[0]]->Fill((*FBarrel)[j].GetEdet(),(*FBarrel)[j].GetRear());
			
			// No punchthrough
			if (((*FBarrel)[j].GetEdet() < 50.) && ((*FBarrel)[j].GetRear() > 500.)) {
			par = 4;
			ForwardBarrel_dEvsLabangle[(*FBarrel)[j].GetID()]->Fill(theta/TMath::Pi()*180.,(*FBarrel)[j].GetRear());	 
			ForwardBarrel_dEvsLabangletot->Fill(theta/TMath::Pi()*180.,(*FBarrel)[j].GetRear()/1000.);	 
		}

			// Particle identification
			if (((*FBarrel)[j].GetEdet() > 50.) && ((*FBarrel)[j].GetRear() > 500.)){
				Ecor = (*FBarrel)[j].GetEdet() + (1. - cos(thetadet))*(*FBarrel)[j].GetRear();
				dEcor = (*FBarrel)[j].GetRear()*cos(thetadet);
				ForwardBarrel_dEvsEcortot->Fill(Ecor/1000.,dEcor/1000.);
				ForwardBarrel_dEvsEcor[(*FBarrel)[j].GetID()]->Fill(Ecor,dEcor);
			 	
				// Protons
				if ((dEcor > for1->Eval(Ecor)) && (dEcor < for2->Eval(Ecor))){
					par = 1;
					ForwardBarrel_dEvsEcorp->Fill(Ecor,dEcor);
					// Energy loss in carbon foil
					Rafter = srim.RMylar[1]->Eval((*FBarrel)[j].GetEdet() + (*FBarrel)[j].GetRear());
					Efor = srim.EMylar[1]->Eval(Rafter + 11.57/cos(thetadet));
					Rafter = srim.RT[1]->Eval(Efor);
					Efor = srim.ET[1]->Eval(Rafter + 5./cos(theta));
					Protons_EtotvsLabangle[(*FBarrel)[j].GetID()]->Fill(theta/TMath::Pi()*180.,Efor);
					Protons_EtotvsLabangletot->Fill(theta/TMath::Pi()*180.,Efor);
					Eex = (1 - m[11]/m[12])*Ein - (1 + m[1]/m[12])*Efor + 2*cos(theta)*sqrt(m[1]*m[11]/(m[12]*m[12]))*sqrt(Efor*Ein) + (m[11] + m[2] - m[12] - m[1]);
					Protons_Eex[(*FBarrel)[j].GetID()]->Fill(Eex);
					Protons_Eextot->Fill(Eex/1000.);
					gamma = (Ein + (m[11] + m[2] - m[12] - m[1]) + Eex - Efor + m[12])/m[12];
					beta = sqrt(1. - 1./(gamma*gamma));
				}

				// Deuterons
				if ((dEcor > for2->Eval(Ecor)) && (dEcor < for3->Eval(Ecor))){
					par = 2;
					ForwardBarrel_dEvsEcord[(*FBarrel)[j].GetID()]->Fill(Ecor,dEcor);
					// Energy loss in carbon foil
					Rafter = srim.RMylar[2]->Eval((*FBarrel)[j].GetEdet() + (*FBarrel)[j].GetRear());
					Efor = srim.EMylar[2]->Eval(Rafter + 11.57/cos(thetadet));
					Rafter = srim.RT[2]->Eval(Efor);
					Efor = srim.ET[2]->Eval(Rafter + 5./cos(theta));
					Deuterons_EtotvsLabangle[(*FBarrel)[j].GetID()]->Fill(theta/TMath::Pi()*180.,Efor);
					Deuterons_EtotvsLabangletot->Fill(theta/TMath::Pi()*180.,Efor);
					Eex = (1 - m[11]/m[11])*Ein - (1 + m[2]/m[11])*Efor + 2*cos(theta)*sqrt(m[2]*m[11]/(m[11]*m[11]))*sqrt(Efor*Ein) + (m[11] + m[2] - m[11] - m[2]);
					Deuterons_Eex[(*FBarrel)[j].GetID()]->Fill(Eex);
					Deuterons_Eextot->Fill(Eex/1000.);
					gamma = (Ein + (m[11] + m[2] - m[11] - m[2]) + Eex - Efor + m[11])/m[11];
					beta = sqrt(1. - 1./(gamma*gamma));
				}

				// Tritons
				if ((dEcor > for3->Eval(Ecor)) && (dEcor < for4->Eval(Ecor))){
					par = 3;
					ForwardBarrel_dEvsEcort->Fill(Ecor,dEcor);
					// Energy loss in carbon foil
					Rafter = srim.RMylar[3]->Eval((*FBarrel)[j].GetEdet() + (*FBarrel)[j].GetRear());
					Efor = srim.EMylar[3]->Eval(Rafter + 11.57/cos(thetadet));
					Rafter = srim.RT[3]->Eval(Efor);
					Efor = srim.ET[3]->Eval(Rafter + 5./cos(theta));
					Tritons_EtotvsLabangle[(*FBarrel)[j].GetID()]->Fill(theta/TMath::Pi()*180.,Efor);
					Tritons_EtotvsLabangletot->Fill(theta/TMath::Pi()*180.,Efor);
					Eex = (1 - m[11]/m[10])*Ein - (1 + m[3]/m[10])*Efor + 2*cos(theta)*sqrt(m[3]*m[11]/(m[10]*m[10]))*sqrt(Efor*Ein) + (m[11] + m[2] - m[10] - m[3]);
					Tritons_Eex[(*FBarrel)[j].GetID()]->Fill(Eex);
					Tritons_Eextot->Fill(Eex/1000.);
					gamma = (Ein + (m[11] + m[2] - m[10] - m[3]) + Eex - Efor + m[10])/m[10];
					beta = sqrt(1. - 1./(gamma*gamma));
				}
			}

			// dE for angle histograms
			for (Int_t k=0;k<20;k++){
				if ((thetadet > (k*.03 + .40)) && (thetadet < ((k+1)*.03 + .40))){
					ForwardBarrel_dEangle[(*FBarrel)[j].GetID()][k]->Fill((*FBarrel)[j].GetRear());
				}
				if ((thetadet > (k*.1 + .4)) && (thetadet < ((k+1)*.1 + .4))){
					ForwardBarrel_dEvsEangle[(*FBarrel)[j].GetID()][k]->Fill((*FBarrel)[j].GetEdet()/1.e3,(*FBarrel)[j].GetRear()/1.e3);
				}
			}
		
			// Energy loss in carbon foil for alphas
			Rafter = srim.RDL[4]->Eval((*FBarrel)[j].GetRear());
			Efor = srim.EDL[4]->Eval(Rafter + .2/cos(thetadet));
			Rafter = srim.RMylar[4]->Eval(Efor);
			Efor = srim.EMylar[4]->Eval(Rafter + 11.57/cos(thetadet));

//			printf("r = %f, E = %f, dE = %f \n",Rafter,Efor,(*FBarrel)[j].GetRear());
			// Efor histograms
			ForwardBarrel_Efor[(*FBarrel)[j].GetID()]->Fill(Efor);
			ForwardBarrel_EforvsDetangle[(*FBarrel)[j].GetID()]->Fill(thetadet,Efor);
			ForwardBarrel_EforvsPos[(*FBarrel)[j].GetID()][(*FBarrel)[j].GetStripNr()[0]]->Fill((*FBarrel)[j].GetStripPos()[0],Efor);
		}
	}
	//end forward barrel
	
	//backward barrel
	for(unsigned int j=0; j<BBarrel->size(); j++){
		if((*BBarrel)[j].GetMult() ==0)
			continue;
		if( (*BBarrel)[j].GetID() < 0 || (*BBarrel)[j].GetID() > 3 ){
			cerr<<"Error in entry "<<i<<": "<<j<<". backward barrel detector id is wrong: "<<(*BBarrel)[j].GetID()<<endl;	
			continue;
		}
		if( (*BBarrel)[j].GetStripNr()[0] < 0 || (*BBarrel)[j].GetStripNr()[0] > 15 ){
			cerr<<"Error in entry "<<i<<": "<<j<<". backward barrel detector strip nr is wrong: "<<(*BBarrel)[j].GetStripNr()[0]<<endl;	
			continue;
		}
	
		par = 5;
		time[(*BBarrel)[j].GetID()+4] = (*BBarrel)[j].GetTime();
      
		//filling histos
		if(all!=0){      
			// Position histograms
			BackwardBarrel_Pos[(*BBarrel)[j].GetID()][(*BBarrel)[j].GetStripNr()[0]]->Fill((*BBarrel)[j].GetStripPos()[0]);
			BackwardBarrel_PosvsStrip[(*BBarrel)[j].GetID()]->Fill((*BBarrel)[j].GetStripNr()[0],(*BBarrel)[j].GetStripPos()[0]);
     	
			// Cylindrical coordinates
			z = 8. + 3.125*((*BBarrel)[j].GetStripNr()[0] + gRandom->Rndm());
			yprime = 29. + 0.*((*BBarrel)[j].GetID()==0) + 2.*((*BBarrel)[j].GetID()==3) - 2.*((*BBarrel)[j].GetID()==1) - 0.*((*BBarrel)[j].GetID()==2);
//			yprime = 29.;
			offset = 0 + .04*((*BBarrel)[j].GetID()==0) + 0.*((*BBarrel)[j].GetID()==1) - .04*((*BBarrel)[j].GetID()==2) - 0.*((*BBarrel)[j].GetID()==3);
//			offset = 0.;
			r = sqrt(yprime*yprime + 50*50*(0.5 + offset - (*BBarrel)[j].GetStripPos()[0])*(0.5 + offset - (*BBarrel)[j].GetStripPos()[0]));
			alpha = atan(((*BBarrel)[j].GetStripPos()[0] - 0.5 - offset)*50/yprime);
			// Spherical coordinates
			theta = atan(r/z) +.005*((*BBarrel)[j].GetID()==0) + 0.*((*BBarrel)[j].GetID()==3) - 0.*((*BBarrel)[j].GetID()==1) - .005*((*BBarrel)[j].GetID()==2);
//			theta = atan(r/z);
			phi = alpha + TMath::Pi()/2.*(*BBarrel)[j].GetID();
			if (phi< 0) phi = phi+2*TMath::Pi();
			// Detector angle
			thetadet = acos(sin(theta)*cos(alpha));
			// Angle histograms
			BackwardBarrel_Theta->Fill(theta);
			BackwardBarrel_Phi->Fill(phi);
			
			// dE histograms
			BackwardBarrel_dE[(*BBarrel)[j].GetID()][(*BBarrel)[j].GetStripNr()[0]]->Fill((*BBarrel)[j].GetRear());
			BackwardBarrel_dEvsPos[(*BBarrel)[j].GetID()]->Fill((*BBarrel)[j].GetStripPos()[0],(*BBarrel)[j].GetRear());
			BackwardBarrel_dEtot[(*BBarrel)[j].GetID()]->Fill((*BBarrel)[j].GetRear());
			BackwardBarrel_dEvsStrip[(*BBarrel)[j].GetID()]->Fill((*BBarrel)[j].GetStripNr()[0],(*BBarrel)[j].GetRear());
			for (Int_t M=0;M<16;M++){
				if (((*BBarrel)[j].GetStripPos()[0] > M*.06667)&&((*BBarrel)[j].GetStripPos()[0] < (M+1)*.06667)){
					BackwardBarrel_dEforPos[(*BBarrel)[j].GetID()][M]->Fill((*BBarrel)[j].GetRear());
				}
			}
			
			// E back histograms
			BackwardBarrel_dEvsEtotal->Fill((*BBarrel)[j].GetEdet(),(*BBarrel)[j].GetRear());
			BackwardBarrel_dEvsEtot[(*BBarrel)[j].GetID()]->Fill((*BBarrel)[j].GetEdet(),(*BBarrel)[j].GetRear());
			BackwardBarrel_dEvsE[(*BBarrel)[j].GetID()][(*BBarrel)[j].GetStripNr()[0]]->Fill((*BBarrel)[j].GetEdet(),(*BBarrel)[j].GetRear());
		
			 // No punchthrough
			 if (((*BBarrel)[j].GetEdet() < 50.) && ((*BBarrel)[j].GetRear() > 500.)) BackwardBarrel_dEvsLabangle[(*BBarrel)[j].GetID()]->Fill(180. - theta/TMath::Pi()*180.,(*BBarrel)[j].GetRear());
			

			// Particle identification
			if (((*BBarrel)[j].GetEdet() > 100.) && ((*BBarrel)[j].GetRear() < 5500.)){				
			 	Ecor = (*BBarrel)[j].GetEdet() + (1. - cos(thetadet))*(*BBarrel)[j].GetRear();
				dEcor = (*BBarrel)[j].GetRear()*cos(thetadet);
				BackwardBarrel_dEvsEcortot->Fill(Ecor,dEcor);
				
				// Protons
				if ((dEcor > back1->Eval(Ecor)) && (dEcor < back2->Eval(Ecor))){
					BackwardBarrel_dEvsEcorp->Fill(Ecor,dEcor);
			 		Protons_EtotvsLabangletot->Fill(-theta/TMath::Pi()*180. + 180.,(*BBarrel)[j].GetRear() + (*BBarrel)[j].GetEdet());
				}
			}
	//		if ((*BBarrel)[j].GetEdet() <= 1000.){
	//		 		Protons_EtotvsLabangle->Fill(-theta/TMath::Pi()*180. + 180.,(*BBarrel)[j].GetRear());
	//		}

			// dE for angle histograms
			for (Int_t k=0;k<6;k++){
				if ((thetadet > (k*.1 + .4)) && (thetadet < ((k+1)*.1 + .4))){
					BackwardBarrel_dEvsEangle[(*BBarrel)[j].GetID()][k]->Fill((*BBarrel)[j].GetEdet(),(*BBarrel)[j].GetRear());
				}
			}
		}
	}
    //end backward barrel
	
	//forward CD
    if(all==1){
		for(unsigned int j=0; j<FCD->size(); j++){
			if( (*FCD)[j].GetID() < 0 || (*FCD)[j].GetID() >= 4 ){
				cerr<<"Error in entry "<<i<<": "<<j<<". backward barrel detector id is wrong: "<<(*FCD)[j].GetID()<<endl;	
				continue;
			}
			par = 6;	
			time[(*FCD)[j].GetID()+8] = (*FCD)[j].GetTime();

			//filling histos
			// Ring dependent histograms
			for( unsigned int k=0;k<(*FCD)[j].GetRingNr().size();k++){
				if(verbose) cout << "ring nr " << (*FCD)[j].GetRingNr()[k] << " with energy " << (*FCD)[j].GetRingEnergy()[k] << endl;
	  
				// dE
				CD_dEvsRing[(*FCD)[j].GetID()]->Fill( (*FCD)[j].GetRingNr()[k], (*FCD)[j].GetRingEnergy()[k] );
				CD_dEforRing[(*FCD)[j].GetID()][(*FCD)[j].GetRingNr()[k]]->Fill((*FCD)[j].GetRingEnergy()[k]);
	
				// E back histograms
				CD_dEvsEtot[(*FCD)[j].GetID()]->Fill((*FCD)[j].GetEdet(),(*FCD)[j].GetRingEnergy()[k]);
				CD_dEvsERing[(*FCD)[j].GetID()][(*FCD)[j].GetRingNr()[k]]->Fill((*FCD)[j].GetEdet(),(*FCD)[j].GetRingEnergy()[k]);
		
				// angle
				theta = atan((9. + (gRandom->Rndm() + (*FCD)[j].GetRingNr()[k])*2.)/64.);
				CD_Theta->Fill(theta);
				CD_dEvsTheta[(*FCD)[j].GetID()]->Fill(theta/TMath::Pi()*180.,(*FCD)[j].GetRingEnergy()[k]/1000.);
				CD_dEvsThetatot->Fill(theta/TMath::Pi()*180.,(*FCD)[j].GetRingEnergy()[k]/1000.);
			
				// 11Be
				if ((*FCD)[j].GetRingEnergy()[k]/1000.>CD1->Eval(theta/TMath::Pi()*180.) && (*FCD)[j].GetRingEnergy()[k]/1000.<CD2->Eval(theta/TMath::Pi()*180.)){
					Be11_Theta[(*FCD)[j].GetID()]->Fill(theta/TMath::Pi()*180.);
				}
			}

			// Strip dependent histograms
			for(unsigned int k=0;k<(*FCD)[j].GetStripNr().size();k++){
				if(verbose) cout << "strip nr " << (*FCD)[j].GetStripNr()[k] << " with energy " << (*FCD)[j].GetStripEnergy()[k] << endl;
			
				CD_dEvsStrip[(*FCD)[j].GetID()]->Fill( (*FCD)[j].GetStripNr()[k], (*FCD)[j].GetStripEnergy()[k] );
				CD_dEforStrip[(*FCD)[j].GetID()][(*FCD)[j].GetStripNr()[k]]->Fill((*FCD)[j].GetStripEnergy()[k]);
				CD_dEvsEStrip[(*FCD)[j].GetID()][(*FCD)[j].GetStripNr()[k]]->Fill((*FCD)[j].GetEdet(),(*FCD)[j].GetStripEnergy()[k]);


				// Angle
				phi = (gRandom->Rndm() + (*FCD)[j].GetStripNr()[k])*.07505 - .52535 - ((*FCD)[j].GetID() - 1)*TMath::Pi()/2.;
				if (phi < 0) phi = phi + 2*TMath::Pi();
				CD_Phi->Fill(phi);
	
				// Strip and ring dependent histograms
				for( unsigned int l=0;l<(*FCD)[j].GetRingNr().size();l++){
					CD_StripvsRing[(*FCD)[j].GetID()]->Fill( (*FCD)[j].GetStripNr()[k], (*FCD)[j].GetRingNr()[l] );
					CD_EStripvsERing[(*FCD)[j].GetID()]->Fill( (*FCD)[j].GetStripEnergy()[k], (*FCD)[j].GetRingEnergy()[l] );
					if (TMath::Abs((*FCD)[j].GetRingEnergy()[l] - (*FCD)[j].GetStripEnergy()[k]) < 500){
						CD_dEtot[(*FCD)[j].GetID()]->Fill((*FCD)[j].GetStripEnergy()[k]);
						CD_EtotvsLabangle->Fill(theta/TMath::Pi()*180.,(*FCD)[j].GetStripEnergy()[k]);
					}
				}
			}
		}
    }
    //end cd
	
	//gammas
    if(all!=2){
		for(unsigned int j=0; j<Miniball->size(); j++){
			if((*Miniball)[j].GetCluID()<0 || (*Miniball)[j].GetCluID()>7){
				cout << "Error in entry "<<i<<": "<<j<<". Germanium detector Cluster id is wrong: "<<(*Miniball)[j].GetCluID()<<endl;	
				continue;
			}

			// Fill histograms
			for(unsigned int k=0; k<(*Miniball)[j].GetCrystal().size();k++){  
				if( (*Miniball)[j].GetCrystal()[k].GetCryID() <0 || (*Miniball)[j].GetCrystal()[k].GetCryID() >2){
					cout << "Error in entry "<<i<<": "<<j<<". Germanium detector Crystal id is wrong: "<<(*Miniball)[j].GetCrystal()[k].GetCryID()<<endl;	
					continue;
				}
			
				Gamma_Etot->Fill((*Miniball)[j].GetCrystal()[k].GetCore());
				Gamma_ECrystals[(*Miniball)[j].GetCluID()][(*Miniball)[j].GetCrystal()[k].GetCryID()]->Fill((*Miniball)[j].GetCrystal()[k].GetCore());
				
				for (Int_t l=0;l<(*Miniball)[j].GetCrystal()[k].GetSegID().size();l++){
					Gamma_ESegments[(*Miniball)[j].GetCluID()][(*Miniball)[j].GetCrystal()[k].GetCryID()][(*Miniball)[j].GetCrystal()[k].GetSegID()[l]]->Fill((*Miniball)[j].GetCrystal()[k].GetCore());
					if (par==1) Gamma_Protons_ESegments[(*Miniball)[j].GetCluID()][(*Miniball)[j].GetCrystal()[k].GetCryID()][(*Miniball)[j].GetCrystal()[k].GetSegID()[l]]->Fill((*Miniball)[j].GetCrystal()[k].GetCore());
					if (par==2) Gamma_Deuterons_ESegments[(*Miniball)[j].GetCluID()][(*Miniball)[j].GetCrystal()[k].GetCryID()][(*Miniball)[j].GetCrystal()[k].GetSegID()[l]]->Fill((*Miniball)[j].GetCrystal()[k].GetCore());
		
					EGamDS = gamma*(*Miniball)[j].GetCrystal()[k].GetCore()*(1. - beta*cos(MBTheta[(*Miniball)[j].GetCluID()][(*Miniball)[j].GetCrystal()[k].GetCryID()][(*Miniball)[j].GetCrystal()[k].GetSegID()[l]]));
					Gamma_EDoppler->Fill(EGamDS);
				
					if (par==1){ 
						Gamma_EProtonDS->Fill(EGamDS);
						Gamma_EexvsEGammaProton->Fill(Eex/1000.,EGamDS);
					}
					if (par==2){ 
						Gamma_EDeuteronDS->Fill(EGamDS);
						Gamma_EexvsEGammaDeuteron->Fill(Eex/1000.,EGamDS);
					}
					if (par==3){
						Gamma_ETritonDS->Fill(EGamDS);
						Gamma_EexvsEGammaTriton->Fill(Eex/1000.,EGamDS);
					}
					if (par==4) Gamma_ENPTDS->Fill(EGamDS);

				}
					if (par==1) Gamma_EProton->Fill((*Miniball)[j].GetCrystal()[k].GetCore());
					if (par==2) Gamma_EDeuteron->Fill((*Miniball)[j].GetCrystal()[k].GetCore());
					if (par==3) Gamma_ETriton->Fill((*Miniball)[j].GetCrystal()[k].GetCore());
					if (par==4) Gamma_ENPT->Fill((*Miniball)[j].GetCrystal()[k].GetCore());
			
				for (Int_t l=0;l<12;l++){
					if (time[l]>0){
						if (par==1) Gamma_EvsdtProton->Fill(25*((*Miniball)[j].GetCrystal()[k].GetTime() - time[l]),(*Miniball)[j].GetCrystal()[k].GetCore()/1000.);
						if (par==2) Gamma_EvsdtDeuteron->Fill(25*((*Miniball)[j].GetCrystal()[k].GetTime() - time[l]),(*Miniball)[j].GetCrystal()[k].GetCore()/1000.);
						if (par==3) Gamma_EvsdtTriton->Fill(25*((*Miniball)[j].GetCrystal()[k].GetTime() - time[l]),(*Miniball)[j].GetCrystal()[k].GetCore()/1000.);
						if (par==4) Gamma_EvsdtNPT->Fill(25*((*Miniball)[j].GetCrystal()[k].GetTime() - time[l]),(*Miniball)[j].GetCrystal()[k].GetCore()/1000.);
						if (par==5) Gamma_EvsdtBack->Fill(25*((*Miniball)[j].GetCrystal()[k].GetTime() - time[l]),(*Miniball)[j].GetCrystal()[k].GetCore()/1000.);
						if (par==6) Gamma_EvsdtCD->Fill(25*((*Miniball)[j].GetCrystal()[k].GetTime() - time[l]),(*Miniball)[j].GetCrystal()[k].GetCore()/1000.);
						Gamma_Evsdttot->Fill(25*((*Miniball)[j].GetCrystal()[k].GetTime() - time[l]),(*Miniball)[j].GetCrystal()[k].GetCore()/1000.);
						if ((*Miniball)[j].GetCrystal()[k].GetCore()>690.&&(*Miniball)[j].GetCrystal()[k].GetCore()<720.) Gamma_dtGe->Fill(25*((*Miniball)[j].GetCrystal()[k].GetTime() - time[l]));
						if ((*Miniball)[j].GetCrystal()[k].GetCore()>2000.&&(*Miniball)[j].GetCrystal()[k].GetCore()<2500.&&(((*Miniball)[j].GetCrystal()[k].GetTime() - time[l])>0)) Gamma_Ecut->Fill((*Miniball)[j].GetCrystal()[k].GetCore());
					}
				}
			}
		}
	}
	//end gammas

    if(i%1000 == 0){
      cout<<setw(5)<<setiosflags(ios::fixed)<<setprecision(1)<<(100.*i)/nentries<<" % done\r"<<flush;
    }
  }

	// Writing histograms
	
	// Single histograms
	// Angles
	ForwardBarrel_Theta->Write("",TObject::kOverwrite);
  	ForwardBarrel_Phi->Write("",TObject::kOverwrite);
  	ForwardBarrel_Theta_det->Write("",TObject::kOverwrite);
  	BackwardBarrel_Theta->Write("",TObject::kOverwrite);
  	BackwardBarrel_Phi->Write("",TObject::kOverwrite);
  	CD_Theta->Write("",TObject::kOverwrite);
  	CD_Phi->Write("",TObject::kOverwrite);
 
	// Energy
	ForwardBarrel_dEvsEtotal->Write("",TObject::kOverwrite);
	BackwardBarrel_dEvsEtotal->Write("",TObject::kOverwrite);
	ForwardBarrel_dEvsEcortot->Write("",TObject::kOverwrite);
	BackwardBarrel_dEvsEcortot->Write("",TObject::kOverwrite);
	ForwardBarrel_dEvsEcorp->Write("",TObject::kOverwrite);
	for (Int_t i=0;i<4;i++) ForwardBarrel_dEvsEcord[i]->Write("",TObject::kOverwrite);
	ForwardBarrel_dEvsEcort->Write("",TObject::kOverwrite);
	BackwardBarrel_dEvsEcorp->Write("",TObject::kOverwrite);

	// Particles
	Protons_EtotvsLabangletot->Write("",TObject::kOverwrite);
	Deuterons_EtotvsLabangletot->Write("",TObject::kOverwrite);
	Tritons_EtotvsLabangletot->Write("",TObject::kOverwrite);
	Protons_Eextot->Write("",TObject::kOverwrite);
	Deuterons_Eextot->Write("",TObject::kOverwrite);
	Tritons_Eextot->Write("",TObject::kOverwrite);
	CD_EtotvsLabangle->Write("",TObject::kOverwrite);
	
	// CD
	// Angle
	CD_dEvsThetatot->Write("",TObject::kOverwrite);
	
	// No Punchthrough
	ForwardBarrel_dEvsLabangletot->Write("",TObject::kOverwrite);
	ForwardBarrel_dEvsLabangle_gammacut->Write("",TObject::kOverwrite);

  	// Pad dependent histograms
	for (Int_t i=0;i<4;i++){
		// Position
		ForwardBarrel_PosvsStrip[i]->Write("",TObject::kOverwrite);
		BackwardBarrel_PosvsStrip[i]->Write("",TObject::kOverwrite);
		
		// dE
		ForwardBarrel_dEtot[i]->Write("",TObject::kOverwrite);
		BackwardBarrel_dEtot[i]->Write("",TObject::kOverwrite);
		ForwardBarrel_dEvsEcor[i]->Write("",TObject::kOverwrite);
		ForwardBarrel_dEvsStrip[i]->Write("",TObject::kOverwrite);
		BackwardBarrel_dEvsStrip[i]->Write("",TObject::kOverwrite);
		BackwardBarrel_dEvsPos[i]->Write("",TObject::kOverwrite);
		ForwardBarrel_dEvsDetangle[i]->Write("",TObject::kOverwrite);

		// Efor
		ForwardBarrel_Efor[i]->Write("",TObject::kOverwrite);
		ForwardBarrel_EforvsDetangle[i]->Write("",TObject::kOverwrite);

		// E backdetector
		ForwardBarrel_dEvsEtot[i]->Write("",TObject::kOverwrite);
		BackwardBarrel_dEvsEtot[i]->Write("",TObject::kOverwrite);

		// No punchthrough
		ForwardBarrel_dEvsLabangle[i]->Write("",TObject::kOverwrite);
		BackwardBarrel_dEvsLabangle[i]->Write("",TObject::kOverwrite);
		
		// Particles
		Protons_EtotvsLabangle[i]->Write("",TObject::kOverwrite);
		Deuterons_EtotvsLabangle[i]->Write("",TObject::kOverwrite);
		Tritons_EtotvsLabangle[i]->Write("",TObject::kOverwrite);
		Protons_Eex[i]->Write("",TObject::kOverwrite);
		Deuterons_Eex[i]->Write("",TObject::kOverwrite);
		Tritons_Eex[i]->Write("",TObject::kOverwrite);

		// strip dependent histograms
		for (Int_t j=0;j<16;j++){
			// Position
			ForwardBarrel_Pos[i][j]->Write("",TObject::kOverwrite);
			BackwardBarrel_Pos[i][j]->Write("",TObject::kOverwrite);
		
			// dE
			ForwardBarrel_dE[i][j]->Write("",TObject::kOverwrite);
			ForwardBarrel_dEvsPos[i][j]->Write("",TObject::kOverwrite);
			for (Int_t k=0;k<16;k++) ForwardBarrel_dEforPos[i][j][k]->Write("",TObject::kOverwrite);
			BackwardBarrel_dE[i][j]->Write("",TObject::kOverwrite);
			BackwardBarrel_dEforPos[i][j]->Write("",TObject::kOverwrite);
			
			// E backdetector
			ForwardBarrel_dEvsE[i][j]->Write("",TObject::kOverwrite);
			BackwardBarrel_dEvsE[i][j]->Write("",TObject::kOverwrite);

			// Efor
			ForwardBarrel_EforvsPos[i][j]->Write("",TObject::kOverwrite);
		}
	
		// Angle dependent histograms
		for (Int_t j=0;j<20;j++){
			// dE
			ForwardBarrel_dEangle[i][j]->Write("",TObject::kOverwrite);
			ForwardBarrel_dEvsEangle[i][j]->Write("",TObject::kOverwrite);
			BackwardBarrel_dEvsEangle[i][j]->Write("",TObject::kOverwrite);
		}
	}
	
	// CD dependent histograms
	for (Int_t i=0;i<4;i++){
		//dE
		CD_dEtot[i]->Write("",TObject::kOverwrite);
		CD_dEvsStrip[i]->Write("",TObject::kOverwrite);
		CD_dEvsRing[i]->Write("",TObject::kOverwrite);
		CD_StripvsRing[i]->Write("",TObject::kOverwrite);
		CD_EStripvsERing[i]->Write("",TObject::kOverwrite);
		CD_dEvsEtot[i]->Write("",TObject::kOverwrite);
		
		// Angle
		CD_dEvsTheta[i]->Write("",TObject::kOverwrite);
		
		// 11Be
		Be11_Theta[i]->Write("",TObject::kOverwrite);

		// Strip/ring dependent histograms
		for (Int_t j=0;j<16;j++){
			CD_dEforStrip[i][j]->Write("",TObject::kOverwrite);
			CD_dEforRing[i][j]->Write("",TObject::kOverwrite);
			CD_dEvsEStrip[i][j]->Write("",TObject::kOverwrite);
			CD_dEvsERing[i][j]->Write("",TObject::kOverwrite);
		}
	}
	
	// Gamma histograms
	Gamma_Etot->Write("",TObject::kOverwrite);
	Gamma_EProton->Write("",TObject::kOverwrite);
	Gamma_EDeuteron->Write("",TObject::kOverwrite);
	Gamma_ETriton->Write("",TObject::kOverwrite);
	Gamma_ENPT->Write("",TObject::kOverwrite);
	Gamma_EProtonDS->Write("",TObject::kOverwrite);
	Gamma_EDeuteronDS->Write("",TObject::kOverwrite);
	Gamma_ETritonDS->Write("",TObject::kOverwrite);
	Gamma_ENPTDS->Write("",TObject::kOverwrite);
	Gamma_EDoppler->Write("",TObject::kOverwrite);
	Gamma_EexvsEGammaProton->Write("",TObject::kOverwrite);
	Gamma_EexvsEGammaDeuteron->Write("",TObject::kOverwrite);
	Gamma_EexvsEGammaTriton->Write("",TObject::kOverwrite);

	// Time
	Gamma_EvsdtProton->Write("",TObject::kOverwrite);
	Gamma_EvsdtDeuteron->Write("",TObject::kOverwrite);
	Gamma_EvsdtTriton->Write("",TObject::kOverwrite);
	Gamma_EvsdtNPT->Write("",TObject::kOverwrite);
	Gamma_EvsdtBack->Write("",TObject::kOverwrite);
	Gamma_EvsdtCD->Write("",TObject::kOverwrite);
	Gamma_Evsdttot->Write("",TObject::kOverwrite);
	Gamma_dtGe->Write("",TObject::kOverwrite);
	Gamma_Ecut->Write("",TObject::kOverwrite);

	// clusters
	for (Int_t i=0;i<8;i++){
		
		// Crystals
		for (Int_t j=0;j<3;j++){
			Gamma_ECrystals[i][j]->Write("",TObject::kOverwrite);
			
			// Segments
			for (Int_t k=0;k<6;k++){
				Gamma_ESegments[i][j][k]->Write("",TObject::kOverwrite);
				Gamma_Protons_ESegments[i][j][k]->Write("",TObject::kOverwrite);
				Gamma_Deuterons_ESegments[i][j][k]->Write("",TObject::kOverwrite);
			}
		}
	}

  outfile->Close();
  delete tr;

  return 0;

}
